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We introduce a generalized Hamiltonian Mean Field Model (gHMF) — XY model with both linear 
and quadratic coupling between spins and explicit Hamiltonian dynamics. In addition to the usual 
paramagnetic and ferromagnetic phases, this model also possesses a nematic phase. The gHMF 
can be solved explicitly using Boltzmann-Gibbs (BG) statistical mechanics, in both canonical and 
microcanonical ensembles. However, when the resulting microcanonical phase diagram is compared 
with the one obtained using molecular dynamics simulations, it is found that the two are very 
different. We will present a dynamical theory which allows us to explicitly calculate the phase 
diagram obtained using molecular dynamics simulations without any adjustable parameters. The 
model illustrates the fundamental role played by dynamics as well the inadequacy of BG statistics 
for systems with long-range forces in the thermodynamic limit. 
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A fundamental concept in statistical mechanics, taught 
in a typical course, is equivalence of ensembles [1]. One 
also learns that mean-field theory becomes exact for sys- 
tems with long-range (LR) interactions 0, However, 
in order to have a well defined thermodynamic limit, in 
this case, special care must be taken. The usual approach 
is to scale the strength of the two-body interaction po- 
tential with the number of particles in the system, N. 
This is the, so-called, Kac prescription — it makes the 
infinitely-long-range two-body interaction infinitesimally 
weak [2|. The thermodynamic limit becomes well de- 
fined, since both the kinetic and the potential contribu- 
tions to the total energy now scale linearly with N, mak- 
ing the energy extensive. Over the last decade, however, 
it has become clear that both the ensemble equivalence 
and the exactness of mean-field theory fail for systems 
with LR interactions The phase-diagrams calcu- 

lated using Boltzmann-Gibbs (BG) statistics in canonical 
and microcanonical ensembles do not always coincide 0] ■ 
Furthermore, molecular dynamics simulations, show that 
isolated LR interacting systems become trapped in quasi- 
stationary states (qSS) the life time of which diverges 
with the number of particles [7- 15|. 

The inapplicability of BG statistics to systems with 
LR forces in thermodynamic limit is a consequence of 
the ergodicity breaking. Scaling of two-body potentials 
with the number of particles — essential for the existence 
of a well defined thermodynamic limit — destroys the 
correlations (collisions) between the particles [17J that 
drive normal short-range interacting systems towards the 
thermodynamic equilibrium. Relaxation to the station- 
ary state of an LR system is, therefore, fundamentally 
different from the collisional (correlational) relaxation 
of normal gases and fluids. Collisionless relaxation re- 
lies on the collective excitations and evaporative cooling 
driven by Landau damping 12J, |18[ . The final stationary 
state reached by a collisionless system is intrinsically non- 
ergodic [lI [l9[ . It does not correspond to the maximum 



of the Boltzmann entropy. To exemplify this dichotomy, 
in this Letter we introduce a new generalized Hamilto- 
nian Mean Field Model (gHMF) — a LR version of the 
model studied in ref. [2q[ . [2lj — which can be solved 
exactly using BG statistical mechanics. We will show 
that the equilibrium phase diagram predicted by the BG 
statistics in the microcanonical ensemble is very differ- 
ent from the one obtained using the molecular dynamics 
(MD) simulations. We will then construct a dynamical 
theory that correctly predicts the location and the order 
of the phase transitions observed in MD simulations. 
The gHMF is described by the Hamiltonian 
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where A € [0,1]. The model can be thought of as either 
XY-spins confined to a line, or as particles restricted to 
move on a circle. The latter interpretation is perhaps 
more convenient when discussing MD simulations with 
equations of motion given by: = dH/dpi and p\ — 
-dH/dOi. 

We define the ferromagnetic and nematic order param- 
eters as mi = jjJ2iLi C0S @i an d m 2 = j? J2iLi cos 2#i, 
respectively. Using the usual statistical mechanics ap- 
proach [B[ , we first calculate the microcanonical entropy 
for the gHMF. 

Within BG statistical mechanics all the thermody- 
namic information is encoded in the phase space volume 
accessible to the system with the total energy E 
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The integral in Eq. ([2]) can be divided into two parts 
kinetic and configurational, 



Q(E,N) = J dKQ kin (K)n conf (E - K) 
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where 
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n conf (E -K)= r H d6 t 6(E -K- U{{0i})), (5) 

J — 7T 

and U is the potential energy, second term in Eq. ([1}. 
Integrating over the momentum degrees of freedom, in 
the thermodynamic limit we obtain 
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The microcanonical entropy per particle is s(s) 



. . 1 1 
s(e) = -ln27r+-+sup 
2 2 K 



-ln2rt+^lnfi con/ (iV(£-«0) 

(7) 

where n = K/N = (E—U)/N = e—u. Since the potential 
energy depends only on mi and 11J2, we define 

r2 m (mi,m 2 )= /^Il^i^dZ 008 ^ _ iVmi) 

x 6(J2cos29i- Nm 2 ), (8) 

which using the Fourier representation of the delta func- 
tion can be written as, 

fl m (m 1 ,m 2 ) = f X oo dx fToodyexp{N [-ixmi 

—iyrri2 + In (J d9 exp(ix cos 9 + iy cos 20)] } . (9) 

The integral can be evaluated using the saddle-point 
method. The extremum corresponds to (x*,y*), which 
must satisfy 
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Defining a = ix* and 6 — iy* and neglecting terms of 
order lower than N, 

lnfl m (mi,m 2 ) = -mia(mi,m 2 ) - m 2 b(m\, m 2 ) 
+ In (J d6* exp [a(mi, m 2 ) cos 6* + b(mi, m 2 ) cos 26*])(12) 

In the thermodynamic limit, we may replace In VL con f(E— 
K) by lnfi m (mi,m 2 ) in Eq. [71 Furthermore, noting 
that k = e — u, where u = (1 — Amf — (1 — A)m|)/2, 
the maximization can be taken with respect to mi,m 2 
instead of k. The entropy per particle then becomes 



s(e) = iln27r+ | 



suPmi.ms [|ln(2 £ - 1 + Am? 
+ (1 — A)m|) — mio(mi, m 2 ) — m 2 6(mi, m 2 ) 
hi (/ d0 exp[a(mi, m 2 ) cos + 6(mi, m 2 ) cos 20]) (13) 
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FIG. 1. Microcanonical phase diagram obtained using BG 
statistics. Solid circles are the two tricritical points. 



with the equilibrium values of the order parameter 
(m\,m 2 ) given by 
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Substituting these expressions into Eqs. (fTU)) and (jlll) . 
we find the equilibrium values of the order parameters 
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where for notational simplicity we have dropped *. In the 
case of a first order phase transition — more than one 
solution of Eqs. (|T6)) and (fT7|) — the equilibrium values of 
mi and m 2 will correspond to the ones that lead to the 
maximum entropy. The resulting microcanonical phase 
diagram is shown in Fig. [T] 

Equation @ requires that the system described by 
the Hamiltonian ^ is ergodic — has equal probabil- 
ity of visiting all possible microstates. To see if this 
is the case, we use Molecular Dynamics (MD) simula- 
tions to study its dynamics. For the gHMF, we are in- 
terested to understand how an ordered (ferromagnetic 
or nematic) state can arise from an originally disor- 
dered homogeneous (paramagnetic) particle distribution 
/o(M = I^ ^ " |0|)Q(Po " The Hamilton's 

equations of motion reduce to a second order differential 
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FIG. 2. The out-of-equilibrium phase diagram of the gHMF. 
The squares and triangles are simulation results for the qSS 
nematic-paramagnetic and para-ferromagnetic phase transi- 
tions, respectively. The shaded area represents the nematic- 
ferromagnetic transition region in which either phase occurs 
with equal probability. To the right of this region, the order is 
ferromagnetic, and to the left, nematic. Black solid lines are 
the theoretical predictions for the transitions. All transitions 
are first order. Insets show the phase space particle distribu- 
tion in different phases. Notice the characteristic core-halo 
structure [l^| both inside nematic and ferromagnetic phases. 
The simulations were performed with N = 10 6 particles for 
the paramagnetic-nematic and paramagnetic-ferromagnetic 
transition, and with N = 10 7 particles to locate the insta- 
bility region between the nematic and ferromagnetic phases. 



equation for 0^, 

Bi = F(0i) (18) 
= -Arai(t) sin0,-(t) - 2(1 - A)m 2 (t) sin 20* (t). 

where F(9) is the force acting on a particle located at 
9, and where we have used the fact that (sin 0(f)) = 
(sin 20(f)) = 0, throughout the dynamical evolution [TBI, 
16j. Comparing the phase diagram obtained using MD 
simulations, we see that it is very different from the pre- 
diction of the microcanonical BG statistical mechanics, 
see Fig. H 

Besides occurring in different regions of the (e, A) 
plane, the phase transitions predicted by the BG statis- 
tics are of the wrong order! While the transitions from 
paramagnetic to ferromagnetic or nematic phases are 
found to be of second order, MD simulations show that 



these transitions are of first order. Furthermore, the sec- 
ond order phase transition line between the nematic and 
the ferromagnetic phase disappears completely and is re- 
placed by a region of instability in which either phase can 
occur with equal probability. 

To understand the results of MD simulations, one must 
forget equilibrium statistical mechanics and return to ki- 
netic theory. In the thermodynamic limit, the dynam- 
ical evolution of the one-particle distribution function 
f(9,p,t) of a system with long-range interactions is gov- 



erned exactly by the Vlasov equation [22J. Vlasov dy- 
namics is collisionless — the relaxation to equilibrium 
comes from Landau damping, a dynamical process in 
which individual particles gain energy from collective os- 
cillations, while the oscillations are damped out. The 
one-particle energy of the gHMF is e = p 2 /2 + 1 — 
Ami cos(0) — (1 — A)m2cos(20). Note that the initial 
particle distribution fo(9,p) has mi = m 2 = 0, so that 
it can be expressed as a function of e. This means that 
fo(9,p) is a stationary solution of the Vlasov equation. 
A phase transition in gHMF, therefore, can occur only 
after a dynamical instability. To explore the non-linear 
stability of the gHMF, we consider a perturbation of the 
initial distribution, such that the maximum momentum 
Po Pm{t) = Po + J2^=o A n(t)cos(n9). We define the 
generalized order parameters as 

m„(f) = (cos(n0)) = J f(9,p,t)coa(n9)dpd9 (19) 

where f(9,p,t) = j^;®^ - \9\)e( Pm (t) - \p\). Note 
that this distribution preserves the phase space density, 
as is required by the Vlasov equation. Performing the 
integration in Eq. (fT9|) . we find that m n (t) = A n (T)/2po. 
Taking two temporal derivatives of m n (t) , we obtain, 



(20) 



-n 2 (p 2 cos(n0)) - n(F(9) sin(n0)), 



where we have used the equation of motion, Eg. (1181) . 
Performing the averages using the distribution function 
/(0,p, f), we obtain the equations of motion for the gen- 
eralized order parameters, 



mi 



12e - 6 - A 



mi = /i(mi,m2,m 3 ,m 4 ) (21) 



m 2 + 2 (12s + A - 7) m 2 = / 2 (mi, m 2 , m 3 , m 4 ) (22) 
m 3 + 27(2e - 1) m 3 = / 3 ("ii, m 2 , m 3 , m 4 ) (23) 
rhi + 48(2e — 1) m 4 = J 4 (mi, m 2 , m 3 , m 4 ) (24) 



where 



fx = mi m 2 1 — 



(A - 1) m 2 m 3 - 

3 (2 e — 1) {m\ + m\ m 3 + m 3 [m 2 (2 + m 2 ) + 2(1 + m 2 ) m 4 ] + 2 mi [m 2 + m 2 , + m 2 + m 2 m 4 + mj } (25) 
/ 2 = A (m 2 — mi m 3 + 2 m 2 m 4 ) — 2 m 2 m 4 — 
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12 (2 e — 1) [nij + m 2 777,4 + 2 mi 7773 (1 + 7772 + 7774) + m\ (1 + 2 7712 + 7774) + 2 7772 + 7774 + TO4 2 )] (26) 

/3=^p [( 2 ~ A ) m2-A 777 4 ]- 

9 (2 e - 1) {mj + 6 to? 7773 + 3 mi [m 2 (2 + 777 2 ) + 2 (1 + m 2 ) m 4 ] + 3 7773 [7773 + 2 (m 2 + 777 2 7774 + 7774)] } (27) 
/ 4 = 2A 77717773 - 4(A - 1)7772 - 

48 (2e - 1) [2mi (1 + m 2 ) 7773 + 777 2 (7772 + 7773) + 2 (m^ + 7773) 7774 + m\ + m? (m 2 + 2m 4 )] (28) 



We have restricted ourselves to the first four generalized 
order parameters, since these are already sufficient to 
understand the phase diagram obtained using MD sim- 
ulations. Note that the right hand sides of Eqs. (|2T1 
- IM|) are non-linear functions, so that the transition 
from paramagnetic-to-ferromagnetic or paramagnetic-to- 
nematic phases is determined by the linear stability of 
these equations. Eqs. (|2"Tj) and (|2"2"j) show that the para- 
magnetic phase becomes unstable to ferromagnetic order- 
ing when 12e — 6 — A < and to nematic ordering when 
12e + A — 7<0. The two stability lines agree perfectly 
with the results of MD simulations, see Fig [21 It is im- 
portant to note that 7773 and 7714 always remain linearly 
stable (recall that e > 0.5 for the initial distribution). 

Linear stability analysis, however, is not sufficient to 
determine the order of the phase transitions for which 
the full non-linear equations must be considered. We first 
note that Eqs. (f2~T1 - l2"4"|) are conservative, they do not ac- 
count for the Landau damping that is responsible for the 
relaxation to equilibrium and formation of the core-halo 
structures .15] , like the ones shown in the insets of Fig [2] 
Phcnomenologically, Landau damping can be included in 
Eqs. (|2T1 - 124|) by introducing terms linear in tt7„. The 
relaxation will then proceed towards the fixed points of 
Eqs. (f2"Tl ~ l2"4")) which can be calculated explicitly. We find 
that when either transition line is crossed, the system 
evolves either to nematic (mi = 0, 7772 7^ 0) or ferromag- 
netic (mi 7^ 0, 7772 + 1 0) fixed points. When crossing the 
paramagnetic-nematic phase transition line, (A < 0.5), 
the order parameter wi\ remains zero, while 7772 jumps by 

V ~ ft ~ °- 459 > independent of A. This theoreti- 
cal prediction is in excellent agreement with the results 
of MD simulation which see a jump in the nematic or- 
der parameter of 0.45, characterizing a strong first-order 
phase transition, see Fig [3] When the paramagnetic- 
ferromagnetic line is crossed (A > 0.5), both 777,1 and 7772 
experience a jump. For A = 0.6, the theory predicts the 
jumps to be 0.5102 and —0.1861, for the ferromagnetic 
and nematic parameters, respectively; while the simula- 
tions find 0.41 and —0.10. For A = 1 the theory predicts 
the respective jumps to be 0.555391 and —0.1129, while 
the simulations find 0.45 and —0.07. It is interesting to 
note that while for the nematic transition the jump in 7772 
is universal — independent of A — for the ferromagnetic 
transition this is not the case. 

What will determine the transition between nematic 




(b) 



0.0 25.0 50.0 75.0 100.0 °^ 2 .0 -1.0 0.0 1.0 2.0 

t p 

FIG. 3. Panel (a) shows the growth and saturation of the 
order parameter 7712 across the paramagnetic-nematic transi- 
tion obtained using MD simulations. The predicted theoreti- 
cal value is 7772 = 0.459, which is in excellent agreement with 
the simulations. In panel (b) the symbols are the momentum 
distribution in the qSS obtained using MD, while the solid 
line depicts the corresponding Maxwell-Boltzmann distribu- 
tion to which the systems should relax in the infinite time 
limit. The parameters are A = 0.2 and u = 0.567. 



and ferromagnetic phases? Deep inside the nematic and 
ferromagnetic phases, Eqs. ([2"T1 - |2"41 possess both sta- 
ble nematic (7771 = 0, 7772 + 1 0) and ferromagnetic fixed 
points (mi + 1 0,7772 + 1 0). Which of these fixed points is 
reached first will depend on the initial condition. Start- 
ing from a paramagnetic distribution fo, in the unstable 
region of the phase diagram, both mi and m 2 will grow 
with time. Eqs. (|21I22I) show that the rate of growth 
of the two order parameters are in general very differ- 
ent, while mi ~ e Alt , where Ai = ^(6 + A - 12e)/2, 
the nematic order parameter grows as 777,2 ~ e A2 *, with 
A2 = \/l4 — 24e — 2A. If the nematic order parameter 
first reaches the value characteristic of the nematic fixed 
point, then nematic order will be established, otherwise 
the phase will be ferromagnetic. Therefore, we expect 
that the nematic-ferromagnetic transition line should be 
given by Ai = A2 (solid line between nematic and ferro- 
magnetic phases in Fig. [2]). This is indeed where the in- 
stability characterizing nematic-to-ferromagnetic region 
is found to be, see Fig. [5] 

We have introduced a generalized Hamiltonian Mean 
Field model. In addition to the usual paramagnetic and 
ferromagnetic phases, this model also possesses a nematic 
phase. We have obtained the phase diagram of the gHMF 
using three different methods: BG statistical mechanics, 
MD simulations, and a new dynamical theory introduced 
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in this paper. The model exemplifies the failure of BG 
statistics to describe isolated systems with LR interac- 
tions, in the thermodynamic limit. This is the first time 
that a complex (multi-phase) out-of-equilibrium phase 
diagram for quasi-stationary states has been calculated 
analytically for a system with LR interactions. 
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FAPERGS, INCT-FCx, and by the US-AFOSR under 
the grant FA9550-12-1-0438. 
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